Imaging a Black Hole using only Closure Quantities

In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 using closure-only imaging. This tutorial is a general introduction to closure-only imaging in Comrade. For an introduction to simultaneous image and instrument modeling, see Stokes I Simultaneous Image and Instrument Modeling

Introduction to Closure Imaging

The EHT is the highest-resolution telescope ever created. Its resolution is equivalent to roughly tracking a hockey puck on the moon when viewing it from the earth. However, the EHT is also a unique interferometer. For one, the data it produces is incredibly sparse. The array is formed from only eight geographic locations around the planet, each with its unique telescope. Additionally, the EHT observes at a much higher frequency than typical interferometers. As a result, it is often difficult to directly provide calibrated data since the source model can be complicated. This implies there can be large instrumental effects often called gains that can corrupt our signal. One way to deal with this is to fit quantities that are independent of gains. These are often called closure quantities. The types of closure quantities are briefly described in Introduction to the VLBI Imaging Problem.

In this tutorial, we will do closure-only modeling of M87 to produce preliminary images of M87.

To get started, we will load Comrade

using Comrade




using Pyehtim

  Activating project at `~/work/Comrade.jl/Comrade.jl/examples`

For reproducibility we use a stable random number genreator

using StableRNGs
rng = StableRNG(123)

StableRNGs.LehmerRNG(state=0x000000000000000000000000000000f7)

Load the Data

To download the data visit https://doi.org/10.25739/g85n-f134 To load the eht-imaging obsdata object we do:

obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "examples", "SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))

Python: <ehtim.obsdata.Obsdata object at 0x7fcf3216b9a0>

Now we do some minor preprocessing:

• Scan average the data since the data have been preprocessed so that the gain phases are coherent.
• Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.

obs = scan_average(obs).add_fractional_noise(0.015)

Python: <ehtim.obsdata.Obsdata object at 0x7fcf32168d60>

Now, we extract our closure quantities from the EHT data set.

dlcamp, dcphase  = extract_table(obs, LogClosureAmplitudes(;snrcut=3), ClosurePhases(;snrcut=3))

(EHTObservation{Float64,Comrade.EHTLogClosureAmplitudeDatum{Float64}, ...}
  source: M87
  mjd: 57849
  frequency: 2.27070703125e11
  bandwidth: 1.856e9
  stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples: 128
, EHTObservation{Float64,Comrade.EHTClosurePhaseDatum{Float64}, ...}
  source: M87
  mjd: 57849
  frequency: 2.27070703125e11
  bandwidth: 1.856e9
  stations: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
  nsamples: 152
)

Build the Model/Posterior

For our model, we will be using an image model that consists of a raster of point sources, convolved with some pulse or kernel to make a ContinuousImage object with it Comrade's. generic image model.

function sky(θ, metadata)
    (;fg, c, σimg) = θ
    (;K, meanpr, grid, cache) = metadata
    # Construct the image model we fix the flux to 0.6 Jy in this case
    cp = meanpr .+ σimg.*c.params
    rast = ((1-fg))*K(to_simplex(CenteredLR(), cp))
    img = IntensityMap(rast, grid)
    m = ContinuousImage(img, cache)
    # Add a large-scale gaussian to deal with the over-resolved mas flux
    g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(fg))
    return m + g
end

sky (generic function with 1 method)

Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of views; typically, 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.

npix = 32
fovxy = μas2rad(150.0)

7.27220521664304e-10

Now, we can feed in the array information to form the cache

grid = imagepixels(fovxy, fovxy, npix, npix)
buffer = IntensityMap(zeros(npix,npix), grid)
cache = create_cache(NFFTAlg(dlcamp), buffer, BSplinePulse{3}())

VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [-4.405690154666661e9 787577.6145833326 … -5.999801315555549e9 -15551.297851562484; -4.523017159111106e9 -1.6838098888888871e6 … 3.059254300444441e9 118294.64453124987]), NFFTPlan with 274 sampling points for an input array of size(32, 32) and an output array of size(274,) with dims 1:2, ComplexF64[0.7014670924816703 - 0.519511453845776im, 0.9999999862106403 - 6.398624051637031e-5im, 0.6947985314160843 - 0.5233374146500747im, 0.9328845732113544 - 0.16336912604935844im, 0.9329421134238738 + 0.16331744240563614im, 0.6948519110447875 + 0.5233083088067162im, 0.7925669922789302 + 0.40478845348863063im, 0.9999999864389528 - 6.519719206334107e-5im, 0.9340555032742773 - 0.1520910238460368im, 0.6896791530194316 - 0.5261728587493705im  …  0.8376666586657548 + 0.1784956664688869im, 0.9196690775195503 + 0.09173305645579018im, 0.9135584686897926 + 0.29199716871629755im, 0.9135610039654034 + 0.2919905620278203im, 0.8899638273254515 + 0.28719490128071723im, 0.9568100109264257 + 0.09814170760184636im, 0.8899622135777404 + 0.287201492222944im, 0.9925668796456906 + 0.0027765725149323534im, 0.8376650833595861 - 0.17848899851834446im, 0.9999999999247222 + 7.335331210317892e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

Now we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior

using VLBIImagePriors, Distributions, DistributionsAD

Since we are using a Gaussian Markov random field prior we need to first specify our mean image. For this work we will use a symmetric Gaussian with a FWHM of 50 μas

fwhmfac = 2*sqrt(2*log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0)./fwhmfac))
imgpr = intensitymap(mpr, grid)

2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   X ∈ 32-element LinRange{Float64,...}
→   Y ∈ 32-element LinRange{Float64,...}
And data, 32×32 NamedDimsArray(::Matrix{Float64}, (:X, :Y)):
                 (-3.52247e-10)  …  (3.29522e-10)  (3.52247e-10)
 (-3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8
 (-3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
 (-3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
 (-2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
 (-2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
 (-2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
    ⋮                            ⋱    ⋮            
  (2.15894e-10)     2.46586e-6        5.12225e-6     2.46586e-6
  (2.38619e-10)     1.51463e-6        3.14629e-6     1.51463e-6
  (2.61345e-10)     8.86092e-7   …    1.84065e-6     8.86092e-7
  (2.84071e-10)     4.93724e-7        1.0256e-6      4.93724e-7
  (3.06796e-10)     2.62014e-7        5.44273e-7     2.62014e-7
  (3.29522e-10)     1.32434e-7        2.751e-7       1.32434e-7
  (3.52247e-10)     6.37537e-8        1.32434e-7     6.37537e-8

Now since we are actually modeling our image on the simplex we need to ensure that our mean image has unit flux

imgpr ./= flux(imgpr)

meanpr = to_real(CenteredLR(), Comrade.baseimage(imgpr))
metadata = (;meanpr,K=CenterImage(imgpr), grid, cache)

(meanpr = [-7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; … ; -6.823167558636962 -6.092113891640146 … -6.092113891640146 -6.823167558636962; -7.55422122563378 -6.823167558636962 … -6.823167558636962 -7.55422122563378], K = VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}([0.9944957386363635 -0.0053267045454545746 … 0.005326704545454565 0.005504261363636378; -0.0053267045454545746 0.9948393969941348 … 0.005160603005865086 0.005326704545454551; … ; 0.005326704545454565 0.005160603005865086 … 0.9948393969941348 -0.005326704545454546; 0.005504261363636378 0.005326704545454551 … -0.005326704545454546 0.9944957386363636], (32, 32)), grid = GriddedKeys{(:X, :Y)}
	X: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
	Y: LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32)
, cache = VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}(VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}(NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 4, 2.0, :kaiser_bessel, AbstractNFFTs.POLYNOMIAL, true, false, true, 0x00000000), [-4.405690154666661e9 787577.6145833326 … -5.999801315555549e9 -15551.297851562484; -4.523017159111106e9 -1.6838098888888871e6 … 3.059254300444441e9 118294.64453124987]), NFFTPlan with 274 sampling points for an input array of size(32, 32) and an output array of size(274,) with dims 1:2, ComplexF64[0.7014670924816703 - 0.519511453845776im, 0.9999999862106403 - 6.398624051637031e-5im, 0.6947985314160843 - 0.5233374146500747im, 0.9328845732113544 - 0.16336912604935844im, 0.9329421134238738 + 0.16331744240563614im, 0.6948519110447875 + 0.5233083088067162im, 0.7925669922789302 + 0.40478845348863063im, 0.9999999864389528 - 6.519719206334107e-5im, 0.9340555032742773 - 0.1520910238460368im, 0.6896791530194316 - 0.5261728587493705im  …  0.8376666586657548 + 0.1784956664688869im, 0.9196690775195503 + 0.09173305645579018im, 0.9135584686897926 + 0.29199716871629755im, 0.9135610039654034 + 0.2919905620278203im, 0.8899638273254515 + 0.28719490128071723im, 0.9568100109264257 + 0.09814170760184636im, 0.8899622135777404 + 0.287201492222944im, 0.9925668796456906 + 0.0027765725149323534im, 0.8376650833595861 - 0.17848899851834446im, 0.9999999999247222 + 7.335331210317892e-6im], BSplinePulse{3}(), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]))

In addition we want a reasonable guess for what the resolution of our image should be. For radio astronomy this is given by roughly the longest baseline in the image. To put this into pixel space we then divide by the pixel size.

beam = beamsize(dlcamp)
rat = (beam/(step(grid.X)))

5.326336637737519

To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. This drastically improves the usual N^3 scaling you get from usual Gaussian Processes.

crcache = MarkovRandomFieldCache(meanpr)

VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}(sparse([1, 2, 32, 33, 993, 1, 2, 3, 34, 994  …  31, 991, 1022, 1023, 1024, 32, 992, 993, 1023, 1024], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 50, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 60, 61, 61, 61, 61, 61, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64, 65, 65, 65, 65, 65, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 84, 84, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86, 86, 86, 86, 87, 87, 87, 87, 87, 88, 88, 88, 88, 88, 89, 89, 89, 89, 89, 90, 90, 90, 90, 90, 91, 91, 91, 91, 91, 92, 92, 92, 92, 92, 93, 93, 93, 93, 93, 94, 94, 94, 94, 94, 95, 95, 95, 95, 95, 96, 96, 96, 96, 96, 97, 97, 97, 97, 97, 98, 98, 98, 98, 98, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 101, 101, 101, 101, 101, 102, 102, 102, 102, 102, 103, 103, 103, 103, 103, 104, 104, 104, 104, 104, 105, 105, 105, 105, 105, 106, 106, 106, 106, 106, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110, 111, 111, 111, 111, 111, 112, 112, 112, 112, 112, 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 116, 116, 116, 116, 116, 117, 117, 117, 117, 117, 118, 118, 118, 118, 118, 119, 119, 119, 119, 119, 120, 120, 120, 120, 120, 121, 121, 121, 121, 121, 122, 122, 122, 122, 122, 123, 123, 123, 123, 123, 124, 124, 124, 124, 124, 125, 125, 125, 125, 125, 126, 126, 126, 126, 126, 127, 127, 127, 127, 127, 128, 128, 128, 128, 128, 129, 129, 129, 129, 129, 130, 130, 130, 130, 130, 131, 131, 131, 131, 131, 132, 132, 132, 132, 132, 133, 133, 133, 133, 133, 134, 134, 134, 134, 134, 135, 135, 135, 135, 135, 136, 136, 136, 136, 136, 137, 137, 137, 137, 137, 138, 138, 138, 138, 138, 139, 139, 139, 139, 139, 140, 140, 140, 140, 140, 141, 141, 141, 141, 141, 142, 142, 142, 142, 142, 143, 143, 143, 143, 143, 144, 144, 144, 144, 144, 145, 145, 145, 145, 145, 146, 146, 146, 146, 146, 147, 147, 147, 147, 147, 148, 148, 148, 148, 148, 149, 149, 149, 149, 149, 150, 150, 150, 150, 150, 151, 151, 151, 151, 151, 152, 152, 152, 152, 152, 153, 153, 153, 153, 153, 154, 154, 154, 154, 154, 155, 155, 155, 155, 155, 156, 156, 156, 156, 156, 157, 157, 157, 157, 157, 158, 158, 158, 158, 158, 159, 159, 159, 159, 159, 160, 160, 160, 160, 160, 161, 161, 161, 161, 161, 162, 162, 162, 162, 162, 163, 163, 163, 163, 163, 164, 164, 164, 164, 164, 165, 165, 165, 165, 165, 166, 166, 166, 166, 166, 167, 167, 167, 167, 167, 168, 168, 168, 168, 168, 169, 169, 169, 169, 169, 170, 170, 170, 170, 170, 171, 171, 171, 171, 171, 172, 172, 172, 172, 172, 173, 173, 173, 173, 173, 174, 174, 174, 174, 174, 175, 175, 175, 175, 175, 176, 176, 176, 176, 176, 177, 177, 177, 177, 177, 178, 178, 178, 178, 178, 179, 179, 179, 179, 179, 180, 180, 180, 180, 180, 181, 181, 181, 181, 181, 182, 182, 182, 182, 182, 183, 183, 183, 183, 183, 184, 184, 184, 184, 184, 185, 185, 185, 185, 185, 186, 186, 186, 186, 186, 187, 187, 187, 187, 187, 188, 188, 188, 188, 188, 189, 189, 189, 189, 189, 190, 190, 190, 190, 190, 191, 191, 191, 191, 191, 192, 192, 192, 192, 192, 193, 193, 193, 193, 193, 194, 194, 194, 194, 194, 195, 195, 195, 195, 195, 196, 196, 196, 196, 196, 197, 197, 197, 197, 197, 198, 198, 198, 198, 198, 199, 199, 199, 199, 199, 200, 200, 200, 200, 200, 201, 201, 201, 201, 201, 202, 202, 202, 202, 202, 203, 203, 203, 203, 203, 204, 204, 204, 204, 204, 205, 205, 205, 205, 205, 206, 206, 206, 206, 206, 207, 207, 207, 207, 207, 208, 208, 208, 208, 208, 209, 209, 209, 209, 209, 210, 210, 210, 210, 210, 211, 211, 211, 211, 211, 212, 212, 212, 212, 212, 213, 213, 213, 213, 213, 214, 214, 214, 214, 214, 215, 215, 215, 215, 215, 216, 216, 216, 216, 216, 217, 217, 217, 217, 217, 218, 218, 218, 218, 218, 219, 219, 219, 219, 219, 220, 220, 220, 220, 220, 221, 221, 221, 221, 221, 222, 222, 222, 222, 222, 223, 223, 223, 223, 223, 224, 224, 224, 224, 224, 225, 225, 225, 225, 225, 226, 226, 226, 226, 226, 227, 227, 227, 227, 227, 228, 228, 228, 228, 228, 229, 229, 229, 229, 229, 230, 230, 230, 230, 230, 231, 231, 231, 231, 231, 232, 232, 232, 232, 232, 233, 233, 233, 233, 233, 234, 234, 234, 234, 234, 235, 235, 235, 235, 235, 236, 236, 236, 236, 236, 237, 237, 237, 237, 237, 238, 238, 238, 238, 238, 239, 239, 239, 239, 239, 240, 240, 240, 240, 240, 241, 241, 241, 241, 241, 242, 242, 242, 242, 242, 243, 243, 243, 243, 243, 244, 244, 244, 244, 244, 245, 245, 245, 245, 245, 246, 246, 246, 246, 246, 247, 247, 247, 247, 247, 248, 248, 248, 248, 248, 249, 249, 249, 249, 249, 250, 250, 250, 250, 250, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 253, 253, 253, 253, 253, 254, 254, 254, 254, 254, 255, 255, 255, 255, 255, 256, 256, 256, 256, 256, 257, 257, 257, 257, 257, 258, 258, 258, 258, 258, 259, 259, 259, 259, 259, 260, 260, 260, 260, 260, 261, 261, 261, 261, 261, 262, 262, 262, 262, 262, 263, 263, 263, 263, 263, 264, 264, 264, 264, 264, 265, 265, 265, 265, 265, 266, 266, 266, 266, 266, 267, 267, 267, 267, 267, 268, 268, 268, 268, 268, 269, 269, 269, 269, 269, 270, 270, 270, 270, 270, 271, 271, 271, 271, 271, 272, 272, 272, 272, 272, 273, 273, 273, 273, 273, 274, 274, 274, 274, 274, 275, 275, 275, 275, 275, 276, 276, 276, 276, 276, 277, 277, 277, 277, 277, 278, 278, 278, 278, 278, 279, 279, 279, 279, 279, 280, 280, 280, 280, 280, 281, 281, 281, 281, 281, 282, 282, 282, 282, 282, 283, 283, 283, 283, 283, 284, 284, 284, 284, 284, 285, 285, 285, 285, 285, 286, 286, 286, 286, 286, 287, 287, 287, 287, 287, 288, 288, 288, 288, 288, 289, 289, 289, 289, 289, 290, 290, 290, 290, 290, 291, 291, 291, 291, 291, 292, 292, 292, 292, 292, 293, 293, 293, 293, 293, 294, 294, 294, 294, 294, 295, 295, 295, 295, 295, 296, 296, 296, 296, 296, 297, 297, 297, 297, 297, 298, 298, 298, 298, 298, 299, 299, 299, 299, 299, 300, 300, 300, 300, 300, 301, 301, 301, 301, 301, 302, 302, 302, 302, 302, 303, 303, 303, 303, 303, 304, 304, 304, 304, 304, 305, 305, 305, 305, 305, 306, 306, 306, 306, 306, 307, 307, 307, 307, 307, 308, 308, 308, 308, 308, 309, 309, 309, 309, 309, 310, 310, 310, 310, 310, 311, 311, 311, 311, 311, 312, 312, 312, 312, 312, 313, 313, 313, 313, 313, 314, 314, 314, 314, 314, 315, 315, 315, 315, 315, 316, 316, 316, 316, 316, 317, 317, 317, 317, 317, 318, 318, 318, 318, 318, 319, 319, 319, 319, 319, 320, 320, 320, 320, 320, 321, 321, 321, 321, 321, 322, 322, 322, 322, 322, 323, 323, 323, 323, 323, 324, 324, 324, 324, 324, 325, 325, 325, 325, 325, 326, 326, 326, 326, 326, 327, 327, 327, 327, 327, 328, 328, 328, 328, 328, 329, 329, 329, 329, 329, 330, 330, 330, 330, 330, 331, 331, 331, 331, 331, 332, 332, 332, 332, 332, 333, 333, 333, 333, 333, 334, 334, 334, 334, 334, 335, 335, 335, 335, 335, 336, 336, 336, 336, 336, 337, 337, 337, 337, 337, 338, 338, 338, 338, 338, 339, 339, 339, 339, 339, 340, 340, 340, 340, 340, 341, 341, 341, 341, 341, 342, 342, 342, 342, 342, 343, 343, 343, 343, 343, 344, 344, 344, 344, 344, 345, 345, 345, 345, 345, 346, 346, 346, 346, 346, 347, 347, 347, 347, 347, 348, 348, 348, 348, 348, 349, 349, 349, 349, 349, 350, 350, 350, 350, 350, 351, 351, 351, 351, 351, 352, 352, 352, 352, 352, 353, 353, 353, 353, 353, 354, 354, 354, 354, 354, 355, 355, 355, 355, 355, 356, 356, 356, 356, 356, 357, 357, 357, 357, 357, 358, 358, 358, 358, 358, 359, 359, 359, 359, 359, 360, 360, 360, 360, 360, 361, 361, 361, 361, 361, 362, 362, 362, 362, 362, 363, 363, 363, 363, 363, 364, 364, 364, 364, 364, 365, 365, 365, 365, 365, 366, 366, 366, 366, 366, 367, 367, 367, 367, 367, 368, 368, 368, 368, 368, 369, 369, 369, 369, 369, 370, 370, 370, 370, 370, 371, 371, 371, 371, 371, 372, 372, 372, 372, 372, 373, 373, 373, 373, 373, 374, 374, 374, 374, 374, 375, 375, 375, 375, 375, 376, 376, 376, 376, 376, 377, 377, 377, 377, 377, 378, 378, 378, 378, 378, 379, 379, 379, 379, 379, 380, 380, 380, 380, 380, 381, 381, 381, 381, 381, 382, 382, 382, 382, 382, 383, 383, 383, 383, 383, 384, 384, 384, 384, 384, 385, 385, 385, 385, 385, 386, 386, 386, 386, 386, 387, 387, 387, 387, 387, 388, 388, 388, 388, 388, 389, 389, 389, 389, 389, 390, 390, 390, 390, 390, 391, 391, 391, 391, 391, 392, 392, 392, 392, 392, 393, 393, 393, 393, 393, 394, 394, 394, 394, 394, 395, 395, 395, 395, 395, 396, 396, 396, 396, 396, 397, 397, 397, 397, 397, 398, 398, 398, 398, 398, 399, 399, 399, 399, 399, 400, 400, 400, 400, 400, 401, 401, 401, 401, 401, 402, 402, 402, 402, 402, 403, 403, 403, 403, 403, 404, 404, 404, 404, 404, 405, 405, 405, 405, 405, 406, 406, 406, 406, 406, 407, 407, 407, 407, 407, 408, 408, 408, 408, 408, 409, 409, 409, 409, 409, 410, 410, 410, 410, 410, 411, 411, 411, 411, 411, 412, 412, 412, 412, 412, 413, 413, 413, 413, 413, 414, 414, 414, 414, 414, 415, 415, 415, 415, 415, 416, 416, 416, 416, 416, 417, 417, 417, 417, 417, 418, 418, 418, 418, 418, 419, 419, 419, 419, 419, 420, 420, 420, 420, 420, 421, 421, 421, 421, 421, 422, 422, 422, 422, 422, 423, 423, 423, 423, 423, 424, 424, 424, 424, 424, 425, 425, 425, 425, 425, 426, 426, 426, 426, 426, 427, 427, 427, 427, 427, 428, 428, 428, 428, 428, 429, 429, 429, 429, 429, 430, 430, 430, 430, 430, 431, 431, 431, 431, 431, 432, 432, 432, 432, 432, 433, 433, 433, 433, 433, 434, 434, 434, 434, 434, 435, 435, 435, 435, 435, 436, 436, 436, 436, 436, 437, 437, 437, 437, 437, 438, 438, 438, 438, 438, 439, 439, 439, 439, 439, 440, 440, 440, 440, 440, 441, 441, 441, 441, 441, 442, 442, 442, 442, 442, 443, 443, 443, 443, 443, 444, 444, 444, 444, 444, 445, 445, 445, 445, 445, 446, 446, 446, 446, 446, 447, 447, 447, 447, 447, 448, 448, 448, 448, 448, 449, 449, 449, 449, 449, 450, 450, 450, 450, 450, 451, 451, 451, 451, 451, 452, 452, 452, 452, 452, 453, 453, 453, 453, 453, 454, 454, 454, 454, 454, 455, 455, 455, 455, 455, 456, 456, 456, 456, 456, 457, 457, 457, 457, 457, 458, 458, 458, 458, 458, 459, 459, 459, 459, 459, 460, 460, 460, 460, 460, 461, 461, 461, 461, 461, 462, 462, 462, 462, 462, 463, 463, 463, 463, 463, 464, 464, 464, 464, 464, 465, 465, 465, 465, 465, 466, 466, 466, 466, 466, 467, 467, 467, 467, 467, 468, 468, 468, 468, 468, 469, 469, 469, 469, 469, 470, 470, 470, 470, 470, 471, 471, 471, 471, 471, 472, 472, 472, 472, 472, 473, 473, 473, 473, 473, 474, 474, 474, 474, 474, 475, 475, 475, 475, 475, 476, 476, 476, 476, 476, 477, 477, 477, 477, 477, 478, 478, 478, 478, 478, 479, 479, 479, 479, 479, 480, 480, 480, 480, 480, 481, 481, 481, 481, 481, 482, 482, 482, 482, 482, 483, 483, 483, 483, 483, 484, 484, 484, 484, 484, 485, 485, 485, 485, 485, 486, 486, 486, 486, 486, 487, 487, 487, 487, 487, 488, 488, 488, 488, 488, 489, 489, 489, 489, 489, 490, 490, 490, 490, 490, 491, 491, 491, 491, 491, 492, 492, 492, 492, 492, 493, 493, 493, 493, 493, 494, 494, 494, 494, 494, 495, 495, 495, 495, 495, 496, 496, 496, 496, 496, 497, 497, 497, 497, 497, 498, 498, 498, 498, 498, 499, 499, 499, 499, 499, 500, 500, 500, 500, 500, 501, 501, 501, 501, 501, 502, 502, 502, 502, 502, 503, 503, 503, 503, 503, 504, 504, 504, 504, 504, 505, 505, 505, 505, 505, 506, 506, 506, 506, 506, 507, 507, 507, 507, 507, 508, 508, 508, 508, 508, 509, 509, 509, 509, 509, 510, 510, 510, 510, 510, 511, 511, 511, 511, 511, 512, 512, 512, 512, 512, 513, 513, 513, 513, 513, 514, 514, 514, 514, 514, 515, 515, 515, 515, 515, 516, 516, 516, 516, 516, 517, 517, 517, 517, 517, 518, 518, 518, 518, 518, 519, 519, 519, 519, 519, 520, 520, 520, 520, 520, 521, 521, 521, 521, 521, 522, 522, 522, 522, 522, 523, 523, 523, 523, 523, 524, 524, 524, 524, 524, 525, 525, 525, 525, 525, 526, 526, 526, 526, 526, 527, 527, 527, 527, 527, 528, 528, 528, 528, 528, 529, 529, 529, 529, 529, 530, 530, 530, 530, 530, 531, 531, 531, 531, 531, 532, 532, 532, 532, 532, 533, 533, 533, 533, 533, 534, 534, 534, 534, 534, 535, 535, 535, 535, 535, 536, 536, 536, 536, 536, 537, 537, 537, 537, 537, 538, 538, 538, 538, 538, 539, 539, 539, 539, 539, 540, 540, 540, 540, 540, 541, 541, 541, 541, 541, 542, 542, 542, 542, 542, 543, 543, 543, 543, 543, 544, 544, 544, 544, 544, 545, 545, 545, 545, 545, 546, 546, 546, 546, 546, 547, 547, 547, 547, 547, 548, 548, 548, 548, 548, 549, 549, 549, 549, 549, 550, 550, 550, 550, 550, 551, 551, 551, 551, 551, 552, 552, 552, 552, 552, 553, 553, 553, 553, 553, 554, 554, 554, 554, 554, 555, 555, 555, 555, 555, 556, 556, 556, 556, 556, 557, 557, 557, 557, 557, 558, 558, 558, 558, 558, 559, 559, 559, 559, 559, 560, 560, 560, 560, 560, 561, 561, 561, 561, 561, 562, 562, 562, 562, 562, 563, 563, 563, 563, 563, 564, 564, 564, 564, 564, 565, 565, 565, 565, 565, 566, 566, 566, 566, 566, 567, 567, 567, 567, 567, 568, 568, 568, 568, 568, 569, 569, 569, 569, 569, 570, 570, 570, 570, 570, 571, 571, 571, 571, 571, 572, 572, 572, 572, 572, 573, 573, 573, 573, 573, 574, 574, 574, 574, 574, 575, 575, 575, 575, 575, 576, 576, 576, 576, 576, 577, 577, 577, 577, 577, 578, 578, 578, 578, 578, 579, 579, 579, 579, 579, 580, 580, 580, 580, 580, 581, 581, 581, 581, 581, 582, 582, 582, 582, 582, 583, 583, 583, 583, 583, 584, 584, 584, 584, 584, 585, 585, 585, 585, 585, 586, 586, 586, 586, 586, 587, 587, 587, 587, 587, 588, 588, 588, 588, 588, 589, 589, 589, 589, 589, 590, 590, 590, 590, 590, 591, 591, 591, 591, 591, 592, 592, 592, 592, 592, 593, 593, 593, 593, 593, 594, 594, 594, 594, 594, 595, 595, 595, 595, 595, 596, 596, 596, 596, 596, 597, 597, 597, 597, 597, 598, 598, 598, 598, 598, 599, 599, 599, 599, 599, 600, 600, 600, 600, 600, 601, 601, 601, 601, 601, 602, 602, 602, 602, 602, 603, 603, 603, 603, 603, 604, 604, 604, 604, 604, 605, 605, 605, 605, 605, 606, 606, 606, 606, 606, 607, 607, 607, 607, 607, 608, 608, 608, 608, 608, 609, 609, 609, 609, 609, 610, 610, 610, 610, 610, 611, 611, 611, 611, 611, 612, 612, 612, 612, 612, 613, 613, 613, 613, 613, 614, 614, 614, 614, 614, 615, 615, 615, 615, 615, 616, 616, 616, 616, 616, 617, 617, 617, 617, 617, 618, 618, 618, 618, 618, 619, 619, 619, 619, 619, 620, 620, 620, 620, 620, 621, 621, 621, 621, 621, 622, 622, 622, 622, 622, 623, 623, 623, 623, 623, 624, 624, 624, 624, 624, 625, 625, 625, 625, 625, 626, 626, 626, 626, 626, 627, 627, 627, 627, 627, 628, 628, 628, 628, 628, 629, 629, 629, 629, 629, 630, 630, 630, 630, 630, 631, 631, 631, 631, 631, 632, 632, 632, 632, 632, 633, 633, 633, 633, 633, 634, 634, 634, 634, 634, 635, 635, 635, 635, 635, 636, 636, 636, 636, 636, 637, 637, 637, 637, 637, 638, 638, 638, 638, 638, 639, 639, 639, 639, 639, 640, 640, 640, 640, 640, 641, 641, 641, 641, 641, 642, 642, 642, 642, 642, 643, 643, 643, 643, 643, 644, 644, 644, 644, 644, 645, 645, 645, 645, 645, 646, 646, 646, 646, 646, 647, 647, 647, 647, 647, 648, 648, 648, 648, 648, 649, 649, 649, 649, 649, 650, 650, 650, 650, 650, 651, 651, 651, 651, 651, 652, 652, 652, 652, 652, 653, 653, 653, 653, 653, 654, 654, 654, 654, 654, 655, 655, 655, 655, 655, 656, 656, 656, 656, 656, 657, 657, 657, 657, 657, 658, 658, 658, 658, 658, 659, 659, 659, 659, 659, 660, 660, 660, 660, 660, 661, 661, 661, 661, 661, 662, 662, 662, 662, 662, 663, 663, 663, 663, 663, 664, 664, 664, 664, 664, 665, 665, 665, 665, 665, 666, 666, 666, 666, 666, 667, 667, 667, 667, 667, 668, 668, 668, 668, 668, 669, 669, 669, 669, 669, 670, 670, 670, 670, 670, 671, 671, 671, 671, 671, 672, 672, 672, 672, 672, 673, 673, 673, 673, 673, 674, 674, 674, 674, 674, 675, 675, 675, 675, 675, 676, 676, 676, 676, 676, 677, 677, 677, 677, 677, 678, 678, 678, 678, 678, 679, 679, 679, 679, 679, 680, 680, 680, 680, 680, 681, 681, 681, 681, 681, 682, 682, 682, 682, 682, 683, 683, 683, 683, 683, 684, 684, 684, 684, 684, 685, 685, 685, 685, 685, 686, 686, 686, 686, 686, 687, 687, 687, 687, 687, 688, 688, 688, 688, 688, 689, 689, 689, 689, 689, 690, 690, 690, 690, 690, 691, 691, 691, 691, 691, 692, 692, 692, 692, 692, 693, 693, 693, 693, 693, 694, 694, 694, 694, 694, 695, 695, 695, 695, 695, 696, 696, 696, 696, 696, 697, 697, 697, 697, 697, 698, 698, 698, 698, 698, 699, 699, 699, 699, 699, 700, 700, 700, 700, 700, 701, 701, 701, 701, 701, 702, 702, 702, 702, 702, 703, 703, 703, 703, 703, 704, 704, 704, 704, 704, 705, 705, 705, 705, 705, 706, 706, 706, 706, 706, 707, 707, 707, 707, 707, 708, 708, 708, 708, 708, 709, 709, 709, 709, 709, 710, 710, 710, 710, 710, 711, 711, 711, 711, 711, 712, 712, 712, 712, 712, 713, 713, 713, 713, 713, 714, 714, 714, 714, 714, 715, 715, 715, 715, 715, 716, 716, 716, 716, 716, 717, 717, 717, 717, 717, 718, 718, 718, 718, 718, 719, 719, 719, 719, 719, 720, 720, 720, 720, 720, 721, 721, 721, 721, 721, 722, 722, 722, 722, 722, 723, 723, 723, 723, 723, 724, 724, 724, 724, 724, 725, 725, 725, 725, 725, 726, 726, 726, 726, 726, 727, 727, 727, 727, 727, 728, 728, 728, 728, 728, 729, 729, 729, 729, 729, 730, 730, 730, 730, 730, 731, 731, 731, 731, 731, 732, 732, 732, 732, 732, 733, 733, 733, 733, 733, 734, 734, 734, 734, 734, 735, 735, 735, 735, 735, 736, 736, 736, 736, 736, 737, 737, 737, 737, 737, 738, 738, 738, 738, 738, 739, 739, 739, 739, 739, 740, 740, 740, 740, 740, 741, 741, 741, 741, 741, 742, 742, 742, 742, 742, 743, 743, 743, 743, 743, 744, 744, 744, 744, 744, 745, 745, 745, 745, 745, 746, 746, 746, 746, 746, 747, 747, 747, 747, 747, 748, 748, 748, 748, 748, 749, 749, 749, 749, 749, 750, 750, 750, 750, 750, 751, 751, 751, 751, 751, 752, 752, 752, 752, 752, 753, 753, 753, 753, 753, 754, 754, 754, 754, 754, 755, 755, 755, 755, 755, 756, 756, 756, 756, 756, 757, 757, 757, 757, 757, 758, 758, 758, 758, 758, 759, 759, 759, 759, 759, 760, 760, 760, 760, 760, 761, 761, 761, 761, 761, 762, 762, 762, 762, 762, 763, 763, 763, 763, 763, 764, 764, 764, 764, 764, 765, 765, 765, 765, 765, 766, 766, 766, 766, 766, 767, 767, 767, 767, 767, 768, 768, 768, 768, 768, 769, 769, 769, 769, 769, 770, 770, 770, 770, 770, 771, 771, 771, 771, 771, 772, 772, 772, 772, 772, 773, 773, 773, 773, 773, 774, 774, 774, 774, 774, 775, 775, 775, 775, 775, 776, 776, 776, 776, 776, 777, 777, 777, 777, 777, 778, 778, 778, 778, 778, 779, 779, 779, 779, 779, 780, 780, 780, 780, 780, 781, 781, 781, 781, 781, 782, 782, 782, 782, 782, 783, 783, 783, 783, 783, 784, 784, 784, 784, 784, 785, 785, 785, 785, 785, 786, 786, 786, 786, 786, 787, 787, 787, 787, 787, 788, 788, 788, 788, 788, 789, 789, 789, 789, 789, 790, 790, 790, 790, 790, 791, 791, 791, 791, 791, 792, 792, 792, 792, 792, 793, 793, 793, 793, 793, 794, 794, 794, 794, 794, 795, 795, 795, 795, 795, 796, 796, 796, 796, 796, 797, 797, 797, 797, 797, 798, 798, 798, 798, 798, 799, 799, 799, 799, 799, 800, 800, 800, 800, 800, 801, 801, 801, 801, 801, 802, 802, 802, 802, 802, 803, 803, 803, 803, 803, 804, 804, 804, 804, 804, 805, 805, 805, 805, 805, 806, 806, 806, 806, 806, 807, 807, 807, 807, 807, 808, 808, 808, 808, 808, 809, 809, 809, 809, 809, 810, 810, 810, 810, 810, 811, 811, 811, 811, 811, 812, 812, 812, 812, 812, 813, 813, 813, 813, 813, 814, 814, 814, 814, 814, 815, 815, 815, 815, 815, 816, 816, 816, 816, 816, 817, 817, 817, 817, 817, 818, 818, 818, 818, 818, 819, 819, 819, 819, 819, 820, 820, 820, 820, 820, 821, 821, 821, 821, 821, 822, 822, 822, 822, 822, 823, 823, 823, 823, 823, 824, 824, 824, 824, 824, 825, 825, 825, 825, 825, 826, 826, 826, 826, 826, 827, 827, 827, 827, 827, 828, 828, 828, 828, 828, 829, 829, 829, 829, 829, 830, 830, 830, 830, 830, 831, 831, 831, 831, 831, 832, 832, 832, 832, 832, 833, 833, 833, 833, 833, 834, 834, 834, 834, 834, 835, 835, 835, 835, 835, 836, 836, 836, 836, 836, 837, 837, 837, 837, 837, 838, 838, 838, 838, 838, 839, 839, 839, 839, 839, 840, 840, 840, 840, 840, 841, 841, 841, 841, 841, 842, 842, 842, 842, 842, 843, 843, 843, 843, 843, 844, 844, 844, 844, 844, 845, 845, 845, 845, 845, 846, 846, 846, 846, 846, 847, 847, 847, 847, 847, 848, 848, 848, 848, 848, 849, 849, 849, 849, 849, 850, 850, 850, 850, 850, 851, 851, 851, 851, 851, 852, 852, 852, 852, 852, 853, 853, 853, 853, 853, 854, 854, 854, 854, 854, 855, 855, 855, 855, 855, 856, 856, 856, 856, 856, 857, 857, 857, 857, 857, 858, 858, 858, 858, 858, 859, 859, 859, 859, 859, 860, 860, 860, 860, 860, 861, 861, 861, 861, 861, 862, 862, 862, 862, 862, 863, 863, 863, 863, 863, 864, 864, 864, 864, 864, 865, 865, 865, 865, 865, 866, 866, 866, 866, 866, 867, 867, 867, 867, 867, 868, 868, 868, 868, 868, 869, 869, 869, 869, 869, 870, 870, 870, 870, 870, 871, 871, 871, 871, 871, 872, 872, 872, 872, 872, 873, 873, 873, 873, 873, 874, 874, 874, 874, 874, 875, 875, 875, 875, 875, 876, 876, 876, 876, 876, 877, 877, 877, 877, 877, 878, 878, 878, 878, 878, 879, 879, 879, 879, 879, 880, 880, 880, 880, 880, 881, 881, 881, 881, 881, 882, 882, 882, 882, 882, 883, 883, 883, 883, 883, 884, 884, 884, 884, 884, 885, 885, 885, 885, 885, 886, 886, 886, 886, 886, 887, 887, 887, 887, 887, 888, 888, 888, 888, 888, 889, 889, 889, 889, 889, 890, 890, 890, 890, 890, 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 893, 893, 893, 893, 894, 894, 894, 894, 894, 895, 895, 895, 895, 895, 896, 896, 896, 896, 896, 897, 897, 897, 897, 897, 898, 898, 898, 898, 898, 899, 899, 899, 899, 899, 900, 900, 900, 900, 900, 901, 901, 901, 901, 901, 902, 902, 902, 902, 902, 903, 903, 903, 903, 903, 904, 904, 904, 904, 904, 905, 905, 905, 905, 905, 906, 906, 906, 906, 906, 907, 907, 907, 907, 907, 908, 908, 908, 908, 908, 909, 909, 909, 909, 909, 910, 910, 910, 910, 910, 911, 911, 911, 911, 911, 912, 912, 912, 912, 912, 913, 913, 913, 913, 913, 914, 914, 914, 914, 914, 915, 915, 915, 915, 915, 916, 916, 916, 916, 916, 917, 917, 917, 917, 917, 918, 918, 918, 918, 918, 919, 919, 919, 919, 919, 920, 920, 920, 920, 920, 921, 921, 921, 921, 921, 922, 922, 922, 922, 922, 923, 923, 923, 923, 923, 924, 924, 924, 924, 924, 925, 925, 925, 925, 925, 926, 926, 926, 926, 926, 927, 927, 927, 927, 927, 928, 928, 928, 928, 928, 929, 929, 929, 929, 929, 930, 930, 930, 930, 930, 931, 931, 931, 931, 931, 932, 932, 932, 932, 932, 933, 933, 933, 933, 933, 934, 934, 934, 934, 934, 935, 935, 935, 935, 935, 936, 936, 936, 936, 936, 937, 937, 937, 937, 937, 938, 938, 938, 938, 938, 939, 939, 939, 939, 939, 940, 940, 940, 940, 940, 941, 941, 941, 941, 941, 942, 942, 942, 942, 942, 943, 943, 943, 943, 943, 944, 944, 944, 944, 944, 945, 945, 945, 945, 945, 946, 946, 946, 946, 946, 947, 947, 947, 947, 947, 948, 948, 948, 948, 948, 949, 949, 949, 949, 949, 950, 950, 950, 950, 950, 951, 951, 951, 951, 951, 952, 952, 952, 952, 952, 953, 953, 953, 953, 953, 954, 954, 954, 954, 954, 955, 955, 955, 955, 955, 956, 956, 956, 956, 956, 957, 957, 957, 957, 957, 958, 958, 958, 958, 958, 959, 959, 959, 959, 959, 960, 960, 960, 960, 960, 961, 961, 961, 961, 961, 962, 962, 962, 962, 962, 963, 963, 963, 963, 963, 964, 964, 964, 964, 964, 965, 965, 965, 965, 965, 966, 966, 966, 966, 966, 967, 967, 967, 967, 967, 968, 968, 968, 968, 968, 969, 969, 969, 969, 969, 970, 970, 970, 970, 970, 971, 971, 971, 971, 971, 972, 972, 972, 972, 972, 973, 973, 973, 973, 973, 974, 974, 974, 974, 974, 975, 975, 975, 975, 975, 976, 976, 976, 976, 976, 977, 977, 977, 977, 977, 978, 978, 978, 978, 978, 979, 979, 979, 979, 979, 980, 980, 980, 980, 980, 981, 981, 981, 981, 981, 982, 982, 982, 982, 982, 983, 983, 983, 983, 983, 984, 984, 984, 984, 984, 985, 985, 985, 985, 985, 986, 986, 986, 986, 986, 987, 987, 987, 987, 987, 988, 988, 988, 988, 988, 989, 989, 989, 989, 989, 990, 990, 990, 990, 990, 991, 991, 991, 991, 991, 992, 992, 992, 992, 992, 993, 993, 993, 993, 993, 994, 994, 994, 994, 994, 995, 995, 995, 995, 995, 996, 996, 996, 996, 996, 997, 997, 997, 997, 997, 998, 998, 998, 998, 998, 999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1001, 1001, 1001, 1001, 1001, 1002, 1002, 1002, 1002, 1002, 1003, 1003, 1003, 1003, 1003, 1004, 1004, 1004, 1004, 1004, 1005, 1005, 1005, 1005, 1005, 1006, 1006, 1006, 1006, 1006, 1007, 1007, 1007, 1007, 1007, 1008, 1008, 1008, 1008, 1008, 1009, 1009, 1009, 1009, 1009, 1010, 1010, 1010, 1010, 1010, 1011, 1011, 1011, 1011, 1011, 1012, 1012, 1012, 1012, 1012, 1013, 1013, 1013, 1013, 1013, 1014, 1014, 1014, 1014, 1014, 1015, 1015, 1015, 1015, 1015, 1016, 1016, 1016, 1016, 1016, 1017, 1017, 1017, 1017, 1017, 1018, 1018, 1018, 1018, 1018, 1019, 1019, 1019, 1019, 1019, 1020, 1020, 1020, 1020, 1020, 1021, 1021, 1021, 1021, 1021, 1022, 1022, 1022, 1022, 1022, 1023, 1023, 1023, 1023, 1023, 1024, 1024, 1024, 1024, 1024], [4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0  …  -1.0, -1.0, -1.0, 4.0, -1.0, -1.0, -1.0, -1.0, -1.0, 4.0], 1024, 1024), [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], [0.0 0.03842943919353914 … 0.15224093497742697 0.03842943919353936; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785; … ; 0.15224093497742697 0.1906703741709661 … 0.30448186995485393 0.19067037417096633; 0.03842943919353914 0.07685887838707828 … 0.1906703741709661 0.0768588783870785])

One of the benefits of the Bayesian approach is that we can fit for the hyperparameters of our prior/regularizers unlike traditional RML appraoches. To construct this heirarchical prior we will first make a map that takes in our regularizer hyperparameters and returns the image prior given those hyperparameters.

fmap = let crcache=crcache
    x->GaussMarkovRandomField(x, 1.0, crcache)
end

#1 (generic function with 1 method)

Now we can finally form our image prior. For this we use a heirarchical prior where the correlation length is given by a inverse gamma prior to prevent overfitting. Gaussian Markov random fields are extremly flexible models. To prevent overfitting it is common to use priors that penalize complexity. Therefore, we want to use priors that enforce similarity to our mean image, and prefer smoothness.

cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.01*rat)))

prior = NamedDist(c = cprior, σimg = truncated(Normal(0.0, 0.1); lower = 0.0), fg=Uniform(0.0, 1.0))

lklhd = RadioLikelihood(sky, dlcamp, dcphase;
                        skymeta = metadata)
post = Posterior(lklhd, prior)

Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:meanpr, :K, :grid, :cache), Tuple{Matrix{Float64}, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}}}}, Nothing, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTLogClosureAmplitudeDatum{Float64}, StructArrays.StructVector{Comrade.EHTLogClosureAmplitudeDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :U4, :V4, :T, :F, :quadrangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{NTuple{4, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}, Comrade.EHTObservation{Float64, Comrade.EHTClosurePhaseDatum{Float64}, StructArrays.StructVector{Comrade.EHTClosurePhaseDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :T, :F, :triangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#34#35"{Float64, Base.Fix2{typeof(logclosure_amplitudes), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}, Comrade.ConditionedLikelihood{Comrade.var"#36#37"{Float64, Base.Fix2{typeof(closure_phases), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}}(RadioLikelihood
	Number of data products: 2
, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3410052124381751)
θ: 2.932506494109095
)

)
, Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0), Distributions.Uniform{Float64}(a=0.0, b=1.0))
)
)

Reconstructing the Image

To sample from this posterior, it is convenient to first move from our constrained parameter space to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This is done using the asflat function.

tpost = asflat(post)

Comrade.TransformedPosterior{Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:meanpr, :K, :grid, :cache), Tuple{Matrix{Float64}, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}}}}, Nothing, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTLogClosureAmplitudeDatum{Float64}, StructArrays.StructVector{Comrade.EHTLogClosureAmplitudeDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :U4, :V4, :T, :F, :quadrangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{NTuple{4, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}, Comrade.EHTObservation{Float64, Comrade.EHTClosurePhaseDatum{Float64}, StructArrays.StructVector{Comrade.EHTClosurePhaseDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :T, :F, :triangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#34#35"{Float64, Base.Fix2{typeof(logclosure_amplitudes), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}, Comrade.ConditionedLikelihood{Comrade.var"#36#37"{Float64, Base.Fix2{typeof(closure_phases), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}}, TransformVariables.TransformTuple{NamedTuple{(:c, :σimg, :fg), Tuple{TransformVariables.TransformTuple{NamedTuple{(:params, :hyperparams), Tuple{TransformVariables.ArrayTransformation{TransformVariables.Identity, 2}, TransformVariables.ShiftedExp{true, Float64}}}}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}}(Posterior{RadioLikelihood{Comrade.ModelMetadata{typeof(Main.sky), NamedTuple{(:meanpr, :K, :grid, :cache), Tuple{Matrix{Float64}, VLBIImagePriors.CenterImage{Matrix{Float64}, Tuple{Int64, Int64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}, VLBISkyModels.NUFTCache{VLBISkyModels.ObservedNUFT{NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}, Matrix{Float64}}, NFFT.NFFTPlan{Float64, 2, 1}, Vector{ComplexF64}, BSplinePulse{3}, KeyedArray{Float64, 2, NamedDimsArray{(:X, :Y), Float64, 2, Matrix{Float64}}, GriddedKeys{(:X, :Y), Tuple{LinRange{Float64, Int64}, LinRange{Float64, Int64}}, ComradeBase.NoHeader, Float64}}}}}}, Nothing, Tuple{Comrade.EHTObservation{Float64, Comrade.EHTLogClosureAmplitudeDatum{Float64}, StructArrays.StructVector{Comrade.EHTLogClosureAmplitudeDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :U4, :V4, :T, :F, :quadrangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{NTuple{4, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}, Comrade.EHTObservation{Float64, Comrade.EHTClosurePhaseDatum{Float64}, StructArrays.StructVector{Comrade.EHTClosurePhaseDatum{Float64}, NamedTuple{(:measurement, :error, :U1, :V1, :U2, :V2, :U3, :V3, :T, :F, :triangle), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol, Symbol}}}}, Int64}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64}}, Tuple{Comrade.ConditionedLikelihood{Comrade.var"#34#35"{Float64, Base.Fix2{typeof(logclosure_amplitudes), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}, Comrade.ConditionedLikelihood{Comrade.var"#36#37"{Float64, Base.Fix2{typeof(closure_phases), Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}, VLBILikelihoods.CholeskyFactor{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.CHOLMOD.Factor{Float64}}}, Vector{Float64}}}, Comrade.ClosureConfig{Comrade.EHTObservation{Float64, Comrade.EHTVisibilityDatum{Float64}, StructArrays.StructVector{Comrade.EHTVisibilityDatum{Float64}, NamedTuple{(:measurement, :error, :U, :V, :T, :F, :baseline), Tuple{Vector{ComplexF64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}}}, Int64}, Comrade.EHTArrayConfiguration{Float64, TypedTables.Table{NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Symbol, Vararg{Float64, 8}}}, 1, NamedTuple{(:sites, :X, :Y, :Z, :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset), Tuple{Vector{Symbol}, Vararg{Vector{Float64}, 8}}}}, TypedTables.Table{NamedTuple{(:start, :stop), Tuple{Float64, Float64}}, 1, NamedTuple{(:start, :stop), Tuple{Vector{Float64}, Vector{Float64}}}}, StructArrays.StructVector{Comrade.ArrayBaselineDatum, NamedTuple{(:U, :V, :T, :F, :baseline, :error, :elevation, :parallactic), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Tuple{Symbol, Symbol}}, Vector{Float64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}, StructArrays.StructVector{Tuple{Float64, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}}}, Int64}}, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, NamedTuple{(:U, :V, :T, :F), NTuple{4, Vector{Float64}}}}, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}}(RadioLikelihood
	Number of data products: 2
, VLBIImagePriors.NamedDist{(:c, :σimg, :fg), Tuple{VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Nothing}, Distributions.Uniform{Float64}}}(
dists: (VLBIImagePriors.HierarchicalPrior{Main.var"#1#2"{VLBIImagePriors.MarkovRandomFieldCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, Distributions.InverseGamma{Float64}}(
priormap: #1
hyperprior: Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.3410052124381751)
θ: 2.932506494109095
)

)
, Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0), Distributions.Uniform{Float64}(a=0.0, b=1.0))
)
), TransformVariables.TransformTuple{NamedTuple{(:c, :σimg, :fg), Tuple{TransformVariables.TransformTuple{NamedTuple{(:params, :hyperparams), Tuple{TransformVariables.ArrayTransformation{TransformVariables.Identity, 2}, TransformVariables.ShiftedExp{true, Float64}}}}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}((c = TransformVariables.TransformTuple{NamedTuple{(:params, :hyperparams), Tuple{TransformVariables.ArrayTransformation{TransformVariables.Identity, 2}, TransformVariables.ShiftedExp{true, Float64}}}}((params = TransformVariables.ArrayTransformation{TransformVariables.Identity, 2}(asℝ, (32, 32)), hyperparams = as(Real, 0.0, ∞)), 1025), σimg = as(Real, 0.0, ∞), fg = as(Real, 0.0, 1.0)), 1027))

We can now also find the dimension of our posterior or the number of parameters we will sample.

ndim = dimension(tpost)

1027

Now we optimize using LBFGS

using ComradeOptimization
using OptimizationOptimJL
using Zygote
f = OptimizationFunction(tpost, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing)
sol = solve(prob, LBFGS(); maxiters=5_00);

Before we analyze our solution we first need to transform back to parameter space.

xopt = transform(tpost, sol)

(c = (params = [-0.7540899656319622 -0.7766682440013806 … -0.7071854085061485 -0.7259413626472166; -0.760526674856007 -0.7802434483682645 … -0.7182192539232248 -0.7357493944621932; … ; -0.7337365131588275 -0.7595175962103428 … -0.6845238237594602 -0.7041350663626422; -0.7445430684518516 -0.7718026773915733 … -0.6961306205740406 -0.7175424129210984], hyperparams = 14.351962592319412), σimg = 0.7477102858287898, fg = 0.03828798506047278)

First we will evaluate our fit by plotting the residuals

using Plots
residual(skymodel(post, xopt), dlcamp, ylabel="Log Closure Amplitude Res.")

and now closure phases

residual(skymodel(post, xopt), dcphase, ylabel="|Closure Phase Res.|")

Now let's plot the MAP estimate.

import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), μas2rad(150.0), μas2rad(150.0), 100, 100)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot, figure=(;resolution=(400, 400),))

To sample from the posterior we will use HMC and more specifically the NUTS algorithm. For information about NUTS see Michael Betancourt's notes.

using ComradeAHMC
using Zygote
metric = DiagEuclideanMetric(ndim)
chain, stats = sample(post, AHMC(;metric, autodiff=Val(:Zygote)), 700; nadapts=500, init_params=xopt)

(NamedTuple{(:c, :σimg, :fg), Tuple{NamedTuple{(:params, :hyperparams), Tuple{Matrix{Float64}, Float64}}, Float64, Float64}}[(c = (params = [-0.7571244711534668 -0.7567103984909943 … -0.7298375804461469 -0.7538111961712505; -0.7505033833362925 -0.789954845538718 … -0.7163146677391313 -0.738499178643207; … ; -0.7385669354420504 -0.7438395593476743 … -0.6304939439850971 -0.7527487546228627; -0.6998062695333326 -0.8469013834445304 … -0.6438575063236519 -0.7024666586025518], hyperparams = 15.072264861769103), σimg = 0.7555113264739167, fg = 0.03842227951020442), (c = (params = [-0.8824516255418571 -0.736876067346221 … -0.7332468930708332 -1.4509681932601255; -0.7185172220776816 -0.9722952191658201 … -1.0403709926795686 -1.618891750450487; … ; -0.5178670101309646 -0.7790753078046321 … -1.0524742956054516 -0.3260445889635208; -1.0288472896227507 -0.868405983476251 … -0.45090463318143814 -1.1841808964158635], hyperparams = 26.783041384142503), σimg = 0.7043347857168719, fg = 0.0919273850016019), (c = (params = [-0.35114541273727384 -0.8409315260785245 … -0.46883849617197704 -1.2375759380694082; -0.1279876200353241 -1.0578190736214466 … -1.3070526477569 -0.8286776260225104; … ; -1.031162179833706 -0.8495899872306393 … -0.5231125512288973 -1.3791808165620838; -1.6862199193038596 -1.2012093557496977 … -0.8163673155442718 -1.5315753128520861], hyperparams = 25.075120020700332), σimg = 0.6969200607894503, fg = 0.08316802760551477), (c = (params = [0.020388569507145512 -0.7546289949925001 … -0.6366170966642783 -1.1478257389366524; -0.36557251284950976 -0.8095060621335206 … -0.169725003029228 -0.3801270776373035; … ; -1.0077489786938492 -0.7646791464872799 … -0.7346056096700039 -0.33616803865999756; -0.43601474378566274 -0.6041251487562171 … -1.298216421213642 -1.839558348234813], hyperparams = 30.41348638483805), σimg = 0.7711335643963944, fg = 0.09798039986922978), (c = (params = [0.19283045590458345 -0.09303279364361015 … -0.9465318818069319 -0.7667715566753696; 0.6703557745853428 0.02570113715601685 … -0.7627316581422919 -0.5724177688538694; … ; -0.9311711284836673 0.24472666714344168 … -0.29243460501440305 -1.1128423796334552; -0.5019452720543586 0.3869628700441396 … -1.0696801103325893 0.1045373922836371], hyperparams = 8.722536536842675), σimg = 0.786919679262737, fg = 0.04664716709448697), (c = (params = [-0.40775647280764965 -0.3728822932244638 … -0.3299310582584553 -1.1477591818628081; -0.7114242505152736 -0.19812515403159037 … -1.9764511755377743 -0.30491611467595503; … ; -0.202818425233409 -1.2233421143554708 … -0.14337974810443382 -0.4615137691119683; -0.8609376221741384 -0.7234514664937248 … -1.1089388996801337 -1.4910982735888236], hyperparams = 23.8922728606703), σimg = 0.775199891781469, fg = 0.04133522941876138), (c = (params = [-1.0820527010096177 -1.3843121172143464 … -0.8193024291179373 -0.7577724494558251; -0.7556836640821453 -0.7538765984482716 … -1.1590403739760378 -0.7281142814301905; … ; -1.236341273611869 0.02952039526965243 … 0.7927083837951566 -0.7781206278352709; -1.509137365912179 -0.9412990065854956 … -0.29091427675393633 -0.6885795389376898], hyperparams = 18.26912662351229), σimg = 0.7641884859084618, fg = 0.08366709710020293), (c = (params = [-1.2755509293237028 -0.31212549638438375 … -1.4491842629405547 -1.1932669067384913; -0.4261413864839489 -0.05688542792979219 … -1.8629590418351067 -0.3414166634487817; … ; -0.06533772605505603 -1.2324535181760337 … -0.848845198910955 -0.48171222999604435; -0.688079222683764 -0.9399200767569476 … -0.547887750136851 -0.5654319260717209], hyperparams = 16.43995566364178), σimg = 0.7020098391507795, fg = 0.09737681143230596), (c = (params = [0.05566245796045683 -0.7548706305072712 … -1.1229008560579974 -1.0103736800173633; -0.6992821478492397 -0.6348312433560426 … -0.07532524774229107 -1.0472639311724627; … ; -1.3394478155837664 -0.5789848886615967 … -0.9631493377160604 -1.61806825458895; -0.9903717522995897 -0.3969079649598854 … -1.4638447693588472 -1.8002611833268587], hyperparams = 30.651492663807208), σimg = 0.6373097727398355, fg = 0.11168517038273994), (c = (params = [-2.1474424111138477 -1.151827708677552 … -0.7421608665533052 -1.41698048903395; -1.0521182579673884 -1.5092804536435094 … -1.9594776462411163 -1.3252115406043026; … ; -0.8864230023874755 -1.339462731379921 … -1.3380125820686322 -1.4159388849038232; -1.7750606333715002 -1.7120363542538284 … -0.9279105949158889 -1.63375363484648], hyperparams = 49.41994110276543), σimg = 0.645262619679076, fg = 0.06419479637381856)  …  (c = (params = [0.0370383700709754 0.310486746568723 … 0.512788644164832 0.5407615083501337; 0.32547573522166334 -0.18116638210481478 … 0.8974757300657248 0.47218577513794097; … ; 0.1592530290804117 -0.11017948700649502 … 0.5163904841088406 -0.46841448847488576; 0.3796200674506581 -0.5116243497899073 … -0.5342641936258473 -1.8557928898347367], hyperparams = 22.758742378953514), σimg = 0.7087078526647314, fg = 0.003699600519184573), (c = (params = [0.050759642668851426 -1.1179033987084266 … -0.3821490653839042 -0.6524654026594408; -0.2248944500713851 -0.4527023215474297 … -0.04266385057048777 -0.22584636088511498; … ; -0.44244983643140756 0.11603839050795035 … 0.292936452809905 -0.13611077182972273; -1.1080306691609818 -0.9906008454923397 … 0.2371041734680177 0.6216534362342462], hyperparams = 24.036681790457756), σimg = 0.7576337047563914, fg = 0.013071038546226224), (c = (params = [-0.15869120617676563 -1.287177392073508 … 0.9371629748996216 -0.4758048964568104; -0.026243558242674513 -0.029664058080486458 … 0.9209674832334283 -0.16960011770175612; … ; -0.24208082994819388 -1.1903542273196135 … 0.21353101700690263 -1.1860315068567955; -0.24860722481354644 -0.08117009651483374 … 0.43834394801393334 -0.43164309494156533], hyperparams = 11.12315829011948), σimg = 0.694689978661194, fg = 0.08378280072884894), (c = (params = [0.6214746174062287 1.4144788807499682 … -1.1545949803345856 -0.010434668035389717; 0.46640006850352844 0.48583469891485964 … -0.7708004120632257 0.5820170403841566; … ; -0.28830655871365407 -0.3826819443403383 … -0.1566564092443885 -0.1628218541558017; -0.6013999464297515 -0.4762071858740237 … -0.8383126343155787 0.47970831997922525], hyperparams = 19.974995496602695), σimg = 0.7188108365877599, fg = 0.034774521606407274), (c = (params = [0.3171264050055294 0.1811385934845696 … 0.6927424256706445 0.08138823577845125; 0.6769891716230769 1.122922073626031 … 0.25026057270166246 -0.6240796465421375; … ; -0.41853790466749013 0.1984350453324647 … -0.9987942021331285 -0.4620310423465807; 0.4624435562288772 1.3710378191369499 … 0.10521950012054596 -0.7960351581849411], hyperparams = 7.290544082388237), σimg = 0.6685707925354557, fg = 0.021507958236788364), (c = (params = [0.21578453526424488 0.1442697648107657 … -1.0239567978502015 0.3439416085988636; -0.9148903673612248 -1.1708224462737906 … 0.10790972251925858 0.020892540132089472; … ; 0.262961518586549 0.5937188037736005 … 0.5882630754381032 0.9354650208854693; 0.4865968420710183 -1.1032311075905172 … 0.5035564541316259 0.950587275310601], hyperparams = 32.01582170857316), σimg = 0.7603541641316821, fg = 0.1389301882632772), (c = (params = [-0.8747916898259362 -0.9733400233059899 … 0.7118214588395114 -0.6019422660461261; -0.7011366839225766 -0.03570892143159917 … -0.8122869308876068 -0.9051316514902017; … ; 0.6829714563321854 -0.10729637383497428 … -0.5290050467766684 -0.058792925745961025; -0.41600266551093845 1.1439647233794024 … -0.3379918007684874 -0.5142865301930875], hyperparams = 16.225230441250535), σimg = 0.7559681497541927, fg = 0.08358644814233594), (c = (params = [-0.1389841837430398 0.7191291196611124 … -0.9386436547992784 -0.5030347199030173; -0.49227443803874893 0.3863444227885982 … -0.3453811203271573 -1.0960418824277431; … ; -0.2107014173264751 0.2651947071865976 … -0.109441706251143 -0.3674721069760209; 0.7111227062988815 -0.14996326151217187 … -0.9581217100392309 -0.5431742060034301], hyperparams = 59.09083334049651), σimg = 0.7412249063256738, fg = 0.176492356878437), (c = (params = [-0.5977068425143846 0.4219293685131459 … 1.1268758853159395 1.4480055146975968; 1.2750380715455611 0.9218897180811993 … 1.4778458221332058 1.1720684873381924; … ; 0.6431664325041712 0.8590407538960824 … 0.05313077647981066 0.29852419705488087; 0.03912463289228977 0.43786820076633587 … 0.7374446827043257 1.0263744352430568], hyperparams = 5.724285356439963), σimg = 0.6676367065587646, fg = 0.009820992457976718), (c = (params = [1.8488373499390915 0.8379499861567681 … 0.4332534207178068 0.3672531660234724; 1.092893991346062 1.5695705130606856 … 0.9340408137506263 1.5670667085674659; … ; -0.03328738974420224 0.2746603358638511 … 0.5446950205638114 0.10762538867122129; 0.3259667401607039 0.4863314654156943 … 0.3664535832807074 0.2649646664488064], hyperparams = 109.87589125405684), σimg = 0.7086933014010586, fg = 0.011027411830099719)], NamedTuple{(:n_steps, :is_accept, :acceptance_rate, :log_density, :hamiltonian_energy, :hamiltonian_energy_error, :max_hamiltonian_energy_error, :tree_depth, :numerical_error, :step_size, :nom_step_size, :is_adapt), Tuple{Int64, Bool, Float64, Float64, Float64, Float64, Float64, Int64, Bool, Float64, Float64, Bool}}[(n_steps = 1023, is_accept = 1, acceptance_rate = 0.9999858563310756, log_density = -189.62870660931085, hamiltonian_energy = 654.0878145988478, hamiltonian_energy_error = 2.57244226986586e-5, max_hamiltonian_energy_error = 2.6241445766572724e-5, tree_depth = 10, numerical_error = 0, step_size = 0.0001, nom_step_size = 0.0001, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9927982999827802, log_density = -620.0384569758219, hamiltonian_energy = 705.2594312645065, hamiltonian_energy_error = 0.007587018001231627, max_hamiltonian_energy_error = 0.014090060264038584, tree_depth = 10, numerical_error = 0, step_size = 0.0015754165898318308, nom_step_size = 0.0015754165898318308, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9991589109025582, log_density = -648.2451009965442, hamiltonian_energy = 1111.2800305960604, hamiltonian_energy_error = -0.009908601514098336, max_hamiltonian_energy_error = -0.03276128215702556, tree_depth = 9, numerical_error = 0, step_size = 0.003194732613566258, nom_step_size = 0.003194732613566258, is_adapt = 1), (n_steps = 255, is_accept = 1, acceptance_rate = 0.4629522758431324, log_density = -656.797312956404, hamiltonian_energy = 1167.2022456833274, hamiltonian_energy_error = -0.169182317124978, max_hamiltonian_energy_error = 2.2189945573013574, tree_depth = 8, numerical_error = 0, step_size = 0.007221414918220456, nom_step_size = 0.007221414918220456, is_adapt = 1), (n_steps = 511, is_accept = 1, acceptance_rate = 0.9165301378779485, log_density = -645.4101478255856, hamiltonian_energy = 1143.13431526501, hamiltonian_energy_error = 0.10504664505947403, max_hamiltonian_energy_error = 0.23672147524439424, tree_depth = 9, numerical_error = 0, step_size = 0.003668199639178762, nom_step_size = 0.003668199639178762, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9567003724570384, log_density = -652.7734113847224, hamiltonian_energy = 1204.8487967135616, hamiltonian_energy_error = -0.1720303238980705, max_hamiltonian_energy_error = -0.46911161309526506, tree_depth = 10, numerical_error = 0, step_size = 0.006377223563061462, nom_step_size = 0.006377223563061462, is_adapt = 1), (n_steps = 511, is_accept = 1, acceptance_rate = 0.5731036494575394, log_density = -697.6844118109107, hamiltonian_energy = 1152.722371041529, hamiltonian_energy_error = -0.29565567743134125, max_hamiltonian_energy_error = 4.732662782934767, tree_depth = 9, numerical_error = 0, step_size = 0.012624094243887406, nom_step_size = 0.012624094243887406, is_adapt = 1), (n_steps = 255, is_accept = 1, acceptance_rate = 0.9587288428630847, log_density = -704.8578348382722, hamiltonian_energy = 1207.0779990645892, hamiltonian_energy_error = 0.07759430016085389, max_hamiltonian_energy_error = 0.14585982048265578, tree_depth = 8, numerical_error = 0, step_size = 0.007591647282653008, nom_step_size = 0.007591647282653008, is_adapt = 1), (n_steps = 127, is_accept = 1, acceptance_rate = 0.19043691700361956, log_density = -691.3413936128472, hamiltonian_energy = 1199.4013328702354, hamiltonian_energy_error = -0.27100473006612447, max_hamiltonian_energy_error = 89.77639986548866, tree_depth = 7, numerical_error = 0, step_size = 0.014918295501868739, nom_step_size = 0.014918295501868739, is_adapt = 1), (n_steps = 1023, is_accept = 1, acceptance_rate = 0.9968448363394408, log_density = -687.034219957132, hamiltonian_energy = 1194.0040595519736, hamiltonian_energy_error = -0.0006374327344929043, max_hamiltonian_energy_error = 0.01980466474719833, tree_depth = 10, numerical_error = 0, step_size = 0.0025821351110419913, nom_step_size = 0.0025821351110419913, is_adapt = 1)  …  (n_steps = 127, is_accept = 1, acceptance_rate = 0.9642360565639214, log_density = -690.5010066158114, hamiltonian_energy = 1204.7735732617566, hamiltonian_energy_error = 0.07974272469141397, max_hamiltonian_energy_error = 0.20041417125867156, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.9942959698233165, log_density = -725.686011173639, hamiltonian_energy = 1228.8948162980246, hamiltonian_energy_error = -0.11889032248154763, max_hamiltonian_energy_error = -0.24585451453208407, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.9982415328526593, log_density = -728.4261535497644, hamiltonian_energy = 1289.9212043675388, hamiltonian_energy_error = -0.016117523487082508, max_hamiltonian_energy_error = -0.13245660081861388, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.9828741556099893, log_density = -737.7201747813074, hamiltonian_energy = 1237.3321426418554, hamiltonian_energy_error = 0.046397571997658815, max_hamiltonian_energy_error = 0.09426451333479235, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.9939509042363277, log_density = -698.6768129294352, hamiltonian_energy = 1245.740227096996, hamiltonian_energy_error = 0.021234462323036496, max_hamiltonian_energy_error = -0.13279021500034105, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.6866657457867021, log_density = -698.7991752448951, hamiltonian_energy = 1203.3631276735607, hamiltonian_energy_error = 0.09467493710803865, max_hamiltonian_energy_error = 0.942409564202535, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.9613424715883185, log_density = -700.0486745193842, hamiltonian_energy = 1218.048885672152, hamiltonian_energy_error = 0.1072697890137988, max_hamiltonian_energy_error = -0.33606116941246, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.9841347648654399, log_density = -697.5376338615272, hamiltonian_energy = 1202.2797721427937, hamiltonian_energy_error = 0.02815719484760848, max_hamiltonian_energy_error = -0.46096954451718375, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.9815072131863001, log_density = -710.5995863826923, hamiltonian_energy = 1219.1607981979312, hamiltonian_energy_error = -0.0014030550944426068, max_hamiltonian_energy_error = -0.436116808306906, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0), (n_steps = 127, is_accept = 1, acceptance_rate = 0.8196473159914998, log_density = -693.5321578359236, hamiltonian_energy = 1216.037859968018, hamiltonian_energy_error = 0.4097318342005565, max_hamiltonian_energy_error = 0.5933677243467628, tree_depth = 7, numerical_error = 0, step_size = 0.02495941659785043, nom_step_size = 0.02495941659785043, is_adapt = 0)])

Now that we have our posterior, we can assess which parts of the image are strongly inferred by the data. This is rather unique to Comrade where more traditional imaging algorithms like CLEAN and RML are inherently unable to assess uncertainty in their reconstructions.

To explore our posterior let's first create images from a bunch of draws from the posterior

msamples = skymodel.(Ref(post), chain[501:2:end]);

The mean image is then given by

using StatsBase
imgs = intensitymap.(msamples, μas2rad(150.0), μas2rad(150.0), 128, 128)
mimg = mean(imgs)
simg = std(imgs)
fig = CM.Figure(;resolution=(400, 400));
CM.image(fig[1,1], mimg,
                   axis=(xreversed=true, aspect=1, title="Mean Image"),
                   colormap=:afmhot)
CM.image(fig[1,2], simg./(max.(mimg, 1e-5)),
                   axis=(xreversed=true, aspect=1, title="1/SNR",), colorrange=(0.0, 2.0),
                   colormap=:afmhot)
CM.image(fig[2,1], imgs[1],
                   axis=(xreversed=true, aspect=1,title="Draw 1"),
                   colormap=:afmhot)
CM.image(fig[2,2], imgs[end],
                   axis=(xreversed=true, aspect=1,title="Draw 2"),
                   colormap=:afmhot)
fig

Now let's see whether our residuals look better.

p = plot();
for s in sample(chain[501:end], 10)
    residual!(p, vlbimodel(post, s), dlcamp)
end
ylabel!("Log-Closure Amplitude Res.");
p

p = plot();
for s in sample(chain[501:end], 10)
    residual!(p, vlbimodel(post, s), dcphase)
end
ylabel!("|Closure Phase Res.|");
p

And viola, you have a quick and preliminary image of M87 fitting only closure products. For a publication-level version we would recommend

1. Running the chain longer and multiple times to properly assess things like ESS and R̂ (see Geometric Modeling of EHT Data)
2. Fitting gains. Typically gain amplitudes are good to 10-20% for the EHT not the infinite uncertainty closures implicitly assume
3. Making sure the posterior is unimodal (hint for this example it isn't!). The EHT image posteriors can be pretty complicated, so typically you want to use a sampler that can deal with multi-modal posteriors. Check out the package Pigeons.jl for an in-development package that should easily enable this type of sampling.

This page was generated using Literate.jl.